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Abstract. Maps of the radial magnetic field at a heliocentric distance of ten solar radii 
are used as boundary conditions in the MHD code CRONOS to simulate a 3D inner- 
heliosphcric solar wind emanating from the rotating Sun out to 1 AU. The input data 
for the magnetic field are the result of solar surface flux transport modelling using ob- 
servational data of sunspot groups coupled with a current-sheet source surface model. 
Amongst several advancements, this allows for higher angular resolution than that of com- 
parable observational data from synoptic magnetograms. The required initial conditions 
for the other MHD quantities are obtained following an empirical approach using an in- 
verse relation between flux tube expansion and radial solar wind speed. The computa- 
tions are performed for representative solar minimum and maximum conditions, and the 
corresponding state of the solar wind up to the Earth's orbit is obtained. After a suc- 
cessful comparison of the latter with observational data, they can be used to drive outer- 
heliosphcric models. 



1. Introduction 

Since Parker's explanation of the expansion of the so- 
lar atmosphere as the solar wind [Parker, 1958], its basic, 
long-term average features have been validated by direct in- 
situ as well as indirect remote measurements. At a radial 
distance of a few solar radii, i.e. at the so-called source sur- 
face, the heliospheric magnetic field (HMF) is essentially 
directed radially and 'frozen in' into the radially expanding 
solar wind plasma. Because the magnetic field remains, due 
to an anchoring of the field line footpoints in the solar atmo- 
sphere, influenced by the solar rotation, the field lines form 
three-dimensional Archimedean or so-called Parker spirals 
[Parker, 1958] if, in a first approximation, the Sun is con- 
sidered as a rigid rotator. 

For shorter periods, a more complex magnetic field struc- 
ture was suggested by Fisk [1996], who took into account 
both the motion of magnetic footpoints on the solar surface 
and the existence of a non-vanishing latitudinal field com- 
ponent. He was able to derive analytical expressions which 
generalize the Parker field correspondingly. Subsequent an- 
alytical investigations [Zurbuchen et al, 1997; Kobylinski, 
2001; Schwadron, 2002; Schwadron et al., 2008] have quan- 
tified the complex structure of this so-called Fisk-field and 
its deviations from the Parker field, and numerical simula- 
tions [Lionello et al, 2006] have revealed that the original 
choice of parameters by Fisk [1996] represented an extreme, 
which has recently been confirmed via a study of the conse- 
quences of the Fisk-field for cosmic ray modulation [Sternal 
et al., 2011]. 

Several authors have suggested modifications of both 
fields [Jokipii and Kota, 1989; Smith and Bieber, 1991; 
Schwadron, 2002; Burger and Hitge, 2004; Burger et al., 
2008; Schwadron et al, 2008], which have been quantita- 
tively compared in a study by Scherer et al. [2010]. While 
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very useful for many purposes, these analytical representa- 
tions of the HMF cannot be used to reproduce actual mea- 
surements in sufficient detail at any given time in the so- 
lar activity cycle. For solar activity minima they are often 
too crude to catch small-scale HMF structures and for the 
HMF during maximum solar activity none of these analyti- 
cal representations is valid. The field expressions suggested 
by Zurbuchen et al. [2004] for use during solar maximum 
are a simple modification of a representation for solar mini- 
mum [Zurbuchen et al, 1997] obtained by assuming a rather 
unlikely so-called Fisk angle, see Sternal et al. [2011]. 

While, physically, the HMF originates in the Sun, it is, 
conceptually, customary to distinguish for modelling pur- 
poses between the coronal magnetic field - filling the re- 
gion from the solar surface out to a spherical 'heliobase' 
[Zhao and Hoeksema, 2010] at several (tens of) solar radii 
- and the HMF beyond. This concept is used in many 
modelling attempts, which can be divided in two groups, 
potential field reconstructions and magnetohydrodynamical 
(MHD) models, see Riley et al. [2006]. The latter approach 
is computationally more expensive but can account for more 
physics, direct time dependence and self-consistency. Exam- 
ples for such MHD modelling of the coronal magnetic field 
are Usmanov and Goldstein [2003], who computed a (tilted) 
axisymmetric solution for solar minimum conditions with 
the assumption of a constant polytropic index, Cohen et al. 
[2007] , who extended this to a varying polytropic index and 
solar maximum conditions, or Lionello et al. [2009]; Riley 
et al. [2011], who employed phenomenological heating func- 
tions that are unrelated to the magnetic field direction. In 
more recent MHD models [e.g., Nakamizo et al, 2009; Feng 
et al, 2012] improved heating and momentum source func- 
tions have been used and the two-region set-up has been 
dropped by exploiting advanced numerical grids providing 
high resolution close to the Sun as well as avoiding both 
coordinate-related singularities and extreme cell size differ- 
ences. For the construction of the heating functions these 
models are, however, implicitly using results from poten- 
tial field reconstructions as, e.g., described in Feng et al. 
[2010]. While, in general, MHD treatments should be pre- 
ferred, all models published so far still have severe limita- 
tions like those indicated for the mentioned examples. 

Not only computationally advantageous, but also with re- 
spect to their ability to resolve structures beyond those that 
can be handled by current MHD models [Riley et al, 2006], 
are potential field models. These pure potential field source 



X - 2 



WIENGARTEN ET AL.: INNER-HELIO SPHERIC MHD SIMULATIONS 



surface (PFSS) models, initially developed by Altachuler 
and Newkirk [1969] and Schatten et al. [1969] were signifi- 
cantly refined by recognizing that the observed photospheric 
field should first be corrected for line-of-sight projection and 
then matched to the radial component of the potential field 
[Wang and Sheeley, 1992], by explicitly taking into account 
additional sheet currents [Zhao and Hoeksema, 1995] re- 
sulting in the so-called current sheet source surface (CSSS) 
model, and by connecting such approaches to the solar wind 
expansion [Arge and Pizzo, 2000], which resulted in the so- 
called Wang-Sheeley-Arge (WSA) model, which is also used 
for the two-region set-up [Pizzo et al, 2011]. 

In the new model presented in this paper we follow a 
similar approach, i.e. we consider a two-region model distin- 
guishing between a coronal magnetic field region and that of 
the HMF which are separated by the heliobase as described 
above. In difference and improvement of earlier approaches, 
however, the input data for the magnetic field are the result 
of solar surface flux transport (SFT) modelling [Jtang et al., 
2010] using observational data of sunspot groups coupled 
with the CSSS model by Zhao and Hoeksema [1995]. This 
approach allows for higher resolution than that of compa- 
rable observational data from synoptic magnetograms. The 
heliobase conditions for the other MHD quantities are ob- 
tained following the empirical approach by Detman et al. 
[2006] and Detman et al. [2011], who employed the WSA 
approach. 

Besides the need for high-resolution models of the three- 
dimensional magnetized solar wind for studies of coronal 
mass ejections [see, e.g. the recent review by Kleimann, 
2012], magnetic clouds [see, e.g. Dalakishvih et al, 2011, 
and references therein] or corotating interaction regions [see, 
e.g. Gosling and Pizzo, 1999], there is the need for improved 
boundary conditions at 1 AU for large-scale heliospheric 
models [see, e.g. Pogorelov et al, 2009]. Our new model 
is of particular interest in the context of long-term mod- 
elling of the heliosphere, because with the employed method 
based on sunspot data the heliobase conditions can be re- 
constructed backward in time for more than three centuries 
[see the recent papers by Jiang et al, 2011a, b]. 

2. The Model 

In its basic setup, CRONOS solves the equations of ideal 
MHD in a one-fluid model, this being a justifiable approach 
for the solar wind plasma [e.g. van der Hoist et al, 2005]. 
In its normalized form (see Table 1) the set of equations to 
be solved reads 
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where p is the mass density, v is the velocity of a fluid 
element, B and E describe the electromagnetic field, e 
is the total energy density, and p is the scalar thermal 
pressure, f is the sum of the gravitational force density 
f g — —pGMo/r 2 e r and additional force terms that may en- 
ter. Furthermore, 1 denotes the unit tensor, and the dyadic 
product is used. So far the system is underdetermined, and 
to complete the set of equations we use the closure relations 
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where 7 is the adiabatic exponent. 

While in the isothermal case (7 = 1) the last term in (5) 
is dropped and equation (3) becomes redundant and needs 



not enter the calculations, a temporally varying tempera- 
ture requires using the full set of equations with 7 / 1. We 
follow Pomoell et al. [2011] in setting 7 = 1.05 to accel- 
erate the solar wind. To model shock structures correctly, 
it would be necessary to work with adiabatic exponents as 
actually present in the solar wind, which are higher. To 
accomplish this, further following the approach of Pomoell 
et al. [2011] we would take the so far obtained solution as 
solar wind background input for a new simulation with in- 
troduced shocks. However, the model presented here does 
not involve shock structures. Nevertheless, extending our 
model in this direction can be a subject for future work. 



Table 1. Summary of normalized quantities and normal- 
ization constants. The choice made here takes length in solar 
radii, and number density and magnetic induction according 
to typical values in the solar photosphere. 
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We use spherical coordinates (r,{),ip) with the origin at 
the center of the Sun. Thus, r is the heliocentric radial 
distance, ■& G [0, it] is the polar angle (with the north pole 
corresponding to $ = 0) and ip G [0, 2-zr] is the azimuthal 
angle. The reference point for ip is as follows: For the test 
cases (see Section 3) azimuthal symmetry is used so that 
results are the same for all ip, and it is not necessary to 
define a reference point. The observationally based data in 
Section 4 is given in Carrington longitudes and the position 
of the Earth serves as a reference point. It depends on the 
time at which the data were taken. This will be described 
in Section 4.1 after explaining how the data were obtained. 
Neglecting differential rotation, solar rotation can be de- 
scribed via an angular frequency 



Q = 14.71°/d 



(8) 



according to Snodgrass and Ulrich [1990], which corresponds 
to the sidereal rotation period of 24.47 days. Accordingly, 
the Sun's rotation axis' tilt of about 7.25° relative to the 
axis of the Earth's orbit is also not taken into account. 
Implementing solar rotation into the model can be done in 
both the solar wind plasma rest-frame and the frame co- 
rotating with the Sun. The first choice introduces azimuthal 
components of B and v by applying respective boundary 
conditions while solving the original set of equations (l)-(7). 
The latter choice cannot be implemented by boundary con- 
ditions only, but requires the introduction of ficititious force 
terms that occur in rotating (and thus non-inertial) frames 
of reference, namely the Coriolis force, the centrifugal force, 
and the Euler force: 

tear = —2pd x v — pfl x (ft x r) — p^ x r. (9) 

at 
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In the chosen spherical coordinate system the angular veloc- 
ity is directed along the z axis of the corresponding Carte- 
sian coordinate system, coinciding with the axis of rotation 
of the Sun: fl = Qe z . Solar rotation is assumed constant 
in time (fl 7^ £l(t)) so that the Euler force term vanishes. 
With the Cartesian unit vector transformed to the spherical 
system this finally gives 



final result is 



f 



(2f2 sin $ Vip + rQ. 2 sin 2 1?) e 



+ (2f2 cos i9 Vtp + rO, sin $ cos 1?) e# 
— (2f2(cos$ v$ + sini? v r )) e v ] . 



(10) 



3. Model Validation: Weber & Davis 
Model 

3.1. Analytical formulation 

Mandatory tests have been performed to check the cor- 
rect implementation of our model. First, Parker's solar 
wind model was reconstructed, using a respective isother- 
mal model. It could be shown that implementations in both 
frames of reference produced correct results of the typical 
radial velocity profile, and magnetic field lines bent to spi- 
rals resulting from solar rotation were obtained. A more 
thorough comparison with analytical results, however, can 
be made with the more complex Weber & Davis model [ We- 
ber and Davis, 1967] (hereafter WD): Parker's assumption 
of frozen-in field lines is not valid close to the Sun since 
there the kinetic energy density is smaller than the magnetic 
energy density. Consequently, it can be assumed that the 
plasma flow follows the magnetic field lines instead well in- 
side the Alfvenic radius ta- This co-rotation should weaken 
when approaching ta, and for the limiting case r ^> va 
Parker's model (v^ oc r" 1 ) should be a good approxima- 
tion. 

Assuming an azimuthal component of the flow velocity and 
solving a steady-state model with axial symmetry in the 
equatorial plane, the final equations of WD for the azimuthal 
components of velocity and magnetic induction in the solar 
wind plasma frame of reference read 



v v (r) = Q.WDT- 



Vr/v(r A ) 
L-Ml 



B v (r) = — .B, 



QwDr 1 - {r/r A ) 2 



(11) 
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v(r A ) 1 

where v(j-a) and Ma denote the radial velocity at the 
Alfvenic critical point and the radial Alfven Mach number, 
respectively. Furthermore, B r = Bo/r 2 as follows from the 
solenoidality condition. It is important to note that, here, 
Slwj) is, according to WD, the "angular velocity of the roots 
of the lines of force in the Sun," whose value is not easily 
accessible. It can be shown [Barker and Marlborough, 1982; 
MacGregor and Pizzo, 1983] that SIwd is related to the ob- 
served angular velocity of the photosphere (as in equation 
(8)) via 

QwD = (l + m (13) 

with 

j, _ Vo,r\Bo,ip\ 

VO,tpBo, r 

where the quantities with zero-subscript indicate the respec- 
tive values at the solar photosphere that are easily accessible 
in a simulation starting at r = _Rq. 

This derivation can be adopted for the co-rotating frame of 
reference as well, taking the initial value for the azimuthal 
velocity v v , co {Rq) = 0, and extending the equation of mo- 
tion by the fictitious force terms (in the ecliptic plane) . The 
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Rq B ^ ( 2 B r (r) 
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B r (r) 



-MiQr 1 - -4 



(15) 




(c) 




Figure 1. (a) The azimuthal velocity (rest frame) ex- 
hibits co-rotation up to the Alfvenic critical point and af- 
terwards decreases as 1/r. It matches the analytic ref- 
erence solution except at the Alfvenic critical point that, 
analytically, is a removable singularity, which cannot be 
resolved here, (b) Azimuthal velocity (co-rotating frame) 
where co-rotation cannot be seen as clearly as for the rest 
frame because the scale is ten times coarser, (c) The az- 
imuthal magnetic field matches the analytic solution only 
for r > 47?0 because of the inner radial boundary condi- 
tion (see text). 
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and 



BAr) 



(RQ/r)vo, r Bo, v + v v (r)B r (r) 
v r (r) 



(16) 



3.2. Numerical simulations 

The analytic results shown above can be compared to the 
numerical results obtained from respective simulations in the 
co-rotating and the rest frame. In both cases the simulation 
box is restricted to extend from the Sun to ten solar radii and 
the number of grid points in radial direction is chosen to be 
N r = 400, yielding a cell size of Ar = 9fl©/400 w 0.027?©. 
Since the WD solution is valid for the ecliptic plane only, 
■& is restricted to the interval [0.47T, 0.67r] and a resolution 
of N$ = 41, using an odd number, gives a cell centered at 
i? = 7r/2 and a cell size of A$ « 0.005-7T. The complete <p 
interval can be covered with a single layer of cells due to the 
symmetry. 

For the implementation of boundary conditions, CRONOS 
makes use of so-called ghost cells that are extensions of 
the actual grid at its boundaries. # boundaries are set to 
"reflecting", which basically ensures conservation of both 
mass- and magnetic flux, but influences the adjacent cells. 
This is compensated by a modest resolution so that there is 
no effect on the equatorial cells. Periodic boundary condi- 
tions are used for tp. The outer radial boundary condition 
is set to "outflow" and is self-explanatory. The inner ra- 
dial boundary condition depends on the respective quantity 
and either fixes it at r — Rq or extrapolates inwards into 
the boundaries ghost cells: For the initial number density, 
n(r) = no(i?©/r) 3 is fixed. Radial velocity is initially set 
to v r (r) = C ■ r with an arbitrary constant C = 0.2 so that 
the steady-state equation of continuity is fulfilled. For this 
flux conserving inward extrapolation is applied. The radial 
magnetic field is initialized as B r (r) = _Bo(-R©/r) 2 . This is 
in accordance with the solenoidality condition since latitu- 
dinal components are set to zero ( B$ = 0, v$ = 0) and the 
azimuthal initialization (see below) yields B v ^ B v (tp). 
The azimuthal components of B and v depend on the frame 
of reference. For the rest frame, v v (r) = Q,r, while in the 
corotating frame, v v%co (r) — initially, both being fixed at 
the inner boundary. The azimuthal magnetic field compo- 
nent is actually not initialized (i.e. set to zero), but the 
computed values at the innermost cells during runtime is 
extrapolated into the inner boundaries ghost cells according 
to 

BAr < Rq) = [BARq + Ar)/B r (R Q + Ar)} ■ B r (r). (17) 

This can be done for both reference frames because the mag- 
netic field is invariant under the respective transformation. 



3.3. Comparison 

Figure 1 shows the results for the azimuthal compo- 
nents for the rest frame and the co-rotating frame: In the 
rest frame (a) v v exhibits the expected behavior as for 
Rq < r < 1.5-R© a decreasing co-rotation is found. In the 
region near the Alfvenic critical radius at r « 2.27?© the in- 
fluence of the magnetic field is still evident: there is not yet 
an 1/r behavior as it develops for even larger r in Parker's 
model. The WD solution (dashed line) matches the code's 
solution except near the Alfvenic point, which numerically 
cannot be exactly resolved due to a singularity that, ana- 
lytically, is a removable one. For the co-rotating frame of 
reference (b) the code also reproduces the analytic solution. 
The co-rotating behavior of the plasma near the solar sur- 
face is only slightly recognizable because the y-axis scale is 
ten times coarser. 

The azimuthal magnetic field (c) is also very similar to that 
of the WD solution, and for r > 4J?© they overlap. The dis- 
crepancy for smaller r can be explained by the inner bound- 
ary condition that extrapolates the ratio between azimuthal 



and radial magnetic field component into the boundary. 
This gives higher absolute values in the boundary because 
B r has an 1/r 2 behavior while for B v it is an 1/r behav- 
ior. These higher absolute values are taken into the simu- 
lation box but their influence vanishes as r increases. Con- 
sequently, a refinement of the boundary condition is not 
essential to get correct results for larger r. 
In summary the code reproduces the analytic solutions of 
WD in both frames of reference while boundary conditions 
can lead to small deviations at the simulation box's bound- 
aries if extrapolation is applied. 

4. Observationally Based Inner Boundary 
Condition 

After successfully reproducing analytic results, we now 
use observationally based input data, i.e. magnetic field 
distributions at ten solar radii, which are the result of a 
solar surface flux transport model coupled with an extrapo- 
lation of the heliospheric field that uses observational data 
of sunspot groups [Jiang et al., 2010]. While this is the key 
input to our code, it is, however, necessary to provide the 
code with initial conditions for all MHD variables. For this 
purpose the approach by Detman et al. [2006, 2011] is fol- 
lowed, who use the empirical inverse relation between flux 
tube expansion factor and solar wind speed to determine the 
latter. 

4.1. Radial magnetic field 

As input data for our modelling, maps of the radial mag- 
netic field at ten solar radii will be used. Their complete 
derivation is described in Jiang et al. [2010] and Schiissler 
and Baumann [2006], but since it is the key input to the 
model presented here, in the following it is summarized 
how the input data was generated. The idea is to use 
sunspot group records as input to a solar surface flux trans- 
port (SFT) model yielding the magnetic flux at the solar 
surface, which is then extrapolated to ten solar radii by 
means of the so-called current sheet source surface (CSSS) 
model: The USAF/NOAA sunspot group records date back 
to 1976 and provide the basis of magnetic flux input, be- 
cause sunspots are associated with emerging bipolar mag- 
netic regions (BMRs). This allows for a model of the whole 
photosphere, whereas in observations, the back of the Sun 
cannot be seen. Along with flux cancellation and transport 
by surface flows the SFT model describes the formation of 
the magnetic flux distribution at the solar surface, which 
then is subject to latitudinal differential rotation, merid- 
ional flow, and turbulent diffusion due to granulation and 
supergranulation, so that all the relevant physical processes 
are implemented. Additionally, computation parameters for 
the magnetic flux are calibrated to match observed values. 
In a next step the HMF has to be determined. The 
most widely used approach to achieve this is the poten- 
tial field source surface (PFSS) model [Schatten et al, 1969; 
Altschuler and Newkirk, 1969], which solves Laplace's equa- 
tion between the photosphere and the so-called source sur- 
face at r = Rss at which the field is forced to be radial. 
Ulysses data, however, motivated an extension of the PFSS 
model to explicitly take into account the heliospheric current 
sheet (HCS), yielding the CSSS model [Zhao and Hoeksema, 
1995; Zhao et al., 2002]. As well as adding horizontal volume 
currents, the CSSS model includes the effect of sheet cur- 
rents by introducing another spherical surface, the so-called 
cusp surface at r = R cs = 1.55-R© < R ss = WRq, within 
which the field contains only horizontal volume currents with 
a characteristic length scale of 0.27?©. This corresponds to 
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Figure 2. Maps of the magnetic field strength at 
(bottom row) at R ss = IORq at solar minimum (in 
according to Jiang et al. [2010]. 

a = 0.2 in the model of Bogdan and Low [1986]. Between 
the cusp surface and the source surface the field is config- 
ured by current sheets so that all field lines passing through 
the cusp surface reach the source surface. The field lines are 
assumed to be radial at the source surface. In our case the 
source surface is located at R ss = 10_Rq and the results of 
the coupled SFT/CSSS model will be used as input to the 
model described here. 

The question as to why the photospheric magnetic field dis- 
tributions are not used to run our MHD model, therefore 
starting at the Sun, can be answered as follows: First, an 
MHD code requires high spatial resolution near the Sun, so 
that even on massively parallel architectures it is time con- 
suming. Additionally, such efforts are known to encounter 
difficult boundary condition related problems, because the 
inner radial boundary lies within the Alfvenic critical radius 
[Nakagawa, 1981a, b]. Second, source surface models are 
easy to implement and require modest computer resources. 
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the solar surface (top row) and the source surface 
1987.2, left) and solar maximum (in 2000.5, right) 



The MHD approach as well as the source surface models 
are both simplifications of reality. Even if MHD models 
incorporate more physics, they cannot claim to cover ev- 
erything. Furthermore, comparisons between the MHD and 
PFSS models have been performed by Riley et al. [2006] 
and their "results endorse the PFSS approach, under the 
right conditions and with appropriate caveats [...]." Given 
that they clearly acknowledge the significance of the var- 
ious refinements of the PFSS model (e.g. current sheets, 
non-spherical source surface) , one should consider their find- 
ing as support not only for the usefulness of PFSS models 
but rather of static extrapolation models in general. Here, 
the input data has been derived using the augmented CSSS 
model, as well as a sophisticated SFT model so that it can 
be assumed that for the general survey presented here the 
magnetic field distributions at ten solar radii are a very good 
input. 

One of the main advantages of the CSSS model is that, as 
discussed in Schiissler and Baumann [2006] , it yields an un- 



X - 6 



WIENGARTEN ET AL.: INNER-HELIOSPHERIC MHD SIMULATIONS 



signed field strength at and beyond the source surface which 
is only very weakly latitude dependent. The PFSS model, 
by contrast has a strong latitudinal variability which is in 
conflict with the Ulysses spacecraft observations. This dif- 
ference in the latitudinal distribution directly affects the ex- 
pansion of field lines from the photosphere to the source 
surface, which is important in our modelling of the plasma 
properties of the inner heliosphere. Hence for our purposes 
the CSSS model is preferred to that of the PFSS model. 
In comparison with synoptic magnetograms often used as 
input data to models such as ours, our input data has some 
advantages since it allows for higher angular resolution and 
the whole photosphere can be modelled whereas in observa- 
tions only part of the Sun is visible. 

Even though a whole time series as output of these models 
exists, we will not make use of all of them in this case study, 
but reserve it for future work. Instead, the focus lies on two 
maps at different times in the solar cycle, i.e. one at so- 
lar minimum (1987.2) and one at solar maximum (2000.5). 
The respective radial magnetic field maps are shown in Fig- 
ure 2. All figures are the same (apart from color scaling) 
as in Figure 4 (a)-(b) and (e)-(f) of Jiang et al. [2010]. 
At solar minimum (left panel) the Sun is relatively quiet, 
and, accordingly, the solar surface (top row) displays few 
to no sunspots and thus BMRs, while the overall structure 
is dipole-like. At the source surface (bottom row) the HCS 
is rather flat and the magnetic field reflects the dipole-like 
structure of the coronal magnetic field becoming radial at 
R ss - However, the overall magnetic field strength is rather 
homogeneous. The latter is also true at solar maximum, but 
the HCS shows strong excursions to high latitudes, and even 
additional current sheets may occur. This is due to the large 
number of BMRs at the solar surface that are more promi- 
nent than the dipole structure. 

The maps cover the full longitudinal and latitudinal intervals 
except for the poles (i.e. i? £ [1°, 179°]) and have a resolu- 
tion of 1° x 1°. The latitudinal coordinate was described in 
Section 2. For the longitudinal coordinate Carrington lon- 
gitudes are used that, together with a given time at which 
the data (which is a snapshot) was taken, define a refer- 
ence system. It is, however, instructive to demonstrate how 
the position of the Earth can be found from the given time 
of the snapshot: The Carrington rotation system of refer- 
ence is based on the synodic solar rotation rate as viewed 
from the Earth, giving the time required for one rotation 
Tcarr = 27.2753 days. At the beginning of a Carrington 
rotation the position of the Earth is at 360° with values de- 
creasing to 0° towards the next rotation. The dates at which 
a new rotation commences can be found in respective tables 
(e.g. http : / /alpo-astronomy . org/ solar/rotn_nos .html). 
Taking the example of the dataset from t Bnap = 1987.2, the 
respective Carrington rotation number is 1786 commencing 
at ten = 1987.15727. The difference in time can be used to 
calculate the longitude of the Earth according to 

VEartH = 360° - 360° • ^ ~ tcR (18) 

which puts the position of the Earth at respective values 
V51987.2 = 154° and 952000.5 = 114° (with Carrington rota- 
tion 1964 commencing at 2000.44899). 

4.2. Remaining MHD quantities 

The simulation requires initial values for the full set of the 
MHD variables, while so far only the radial magnetic field 
component is known from the previously described data set. 
One possible approach to derive the remaining quantities 
from the radial magnetic field data at the inner radial grid 
boundary is described in Detman et al. [2006, 2011] and shall 
be briefly summarized: 



As input data a sequence of photospheric magnetic maps 
composed of daily magnetograms is used. A source surface 
current sheet model is used to extrapolate the HMF, which 
is then fed into an MHD code with an inner radial grid 
boundary at R g t = 0.1AU = 21.5-Rq. It is not possible to 
start the MHD simulations at the source surface most of the 
time. Instead, the Alfvenic critical point must be located 
outside the grid, which ensures that the plasma flow speed 
is higher than any characteristic plasma speed (specifically 
the fast mode wave speed), so that there cannot be any flow 
of information back into the boundary. The lower radial grid 
boundary is thus set to 0.1 AU, where this condition is met. 
It is then required to develop an interface to translate the 
input data from the source surface model to the lower radial 
grid boundary, which gives maps of all the necessary MHD 
variables there. This can only be achieved by empirical for- 
mulas, because it has to mimic the main plasma physics 
in the corona, which is not incorporated in source surface 
models. The resulting set of equations contains a number 
of adjustable parameters that can be tuned to give best fits 
to observational data at Earth. This is an ongoing process 
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Figure 3. Location of footpoints of open magnetic field 
lines in the photosphere (black/red crosses) for minimum 
(a) and maximum (b) conditions. In (a) yellow dots indi- 
cate coronal hole boundaries and cyan dots indicate field 
lines with the smallest expansion factors and thus result 
in the highest velocities. In (b) a comparison with coronal 
holes inferred from National Solar Observatory He maps 
(black contours) is superimposed. The indicated polarities 
were also found to be in good agreement with our data. 
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Figure 4. Initial values at the inner radial grid boundary (left panels) and results at the outer radial 
grid boundary for the converged state (right panels) for solar minimum. The shown quantities are (top 
to bottom) radial velocity, mass density, radial- and azimuthal magnetic field. For the nomalizations see 
Table 1. 



and will go hand in hand with advancements in other areas 
to further augment these kinds of models. The following 
set of equations is taken from Detman et al. [2006] , but has 
been adapted: An augmented formula for the radial veloc- 
ity is taken from Detman et al. [2011], which additionally 



uses the footpoint distance dFP to the nearest coronal hole 
boundary because using both footpoint distance and expan- 
sion factor f 3 gives better results than either one alone [Arge 
et al., 2003], which was found to be true for the data used 
here as well. Furthermore, the azimuthal velocity is trans- 
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formed to the co-rotating reference frame. The following set 
of equations is to be understood as to give initial values at 
the inner radial grid boundary at R gb , while initialisitation 
at larger radii will be adressed subsequently: 
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Before describing the parameters in and the implementation 
of these formulas, it should be mentioned that there are dif- 
ferences between this model and Detman et al. [2006, 2011]. 
Most importantly, their extrapolation scheme puts the 
source surface at smaller radial distances, and the extrapo- 
lation method is not the same as the one that generated the 
input data used here. Thus, their parameter tuning might 
not be optimal for our purpose, but repeating this process 
for the data used here would go beyond the scope of this 
work. 

Equation (19) is based on the empirical inverse corre- 
lation between flux tube expansion factor f s (hereafter EF) 
and radial solar wind speed v r observed near Earth reported 
for potential field models as described in Wang and Sheeley 
[1990] and Wang et al. [1997]. Additionally, d FP is the foot- 
point distance to the nearest coronal hole boundary. The 
remaining unknowns are tuning parameters. The formula 
for the expansion factor reads 



-Re 

Rss 



B(Ro,d ,ip ) 
B(R SSl d,<p) 



(27) 



where B(Rq, $o, <Po) and B(R sa , <p) are the magnetic field 
strengths at a point in the photosphere and (following a 
magnetic field line) at the source surface at R S3 . In or- 
der to know how points on the two respective magnetic 
field maps are connected, it is necessary to find a mapping 
if) H» (#o,¥>o)- The resulting footpoint locations are 
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shown in Figure 3 as black/red dots. During solar minimum 
(a) field lines originate from polar regions, but in this case 
there also exist excursions extending towards the equator. 
At solar maximum (b) there is no discernible structure to be 
made out: Field lines originate from coronal holes that can 
be located anywhere on the solar surface. To test whether 
these results are realistic, they can be compared to their 
observational counterparts, which are coronal holes inferred 
from National Solar Observatory He maps. This is shown 
for the case of solar maximum in Figure 3: Model data 
is shown in red dots, observational data in black contours. 
There is some discrepancy, because the model is based on 
very different data (sunspot areas) and is a snapshot, while 
the observed synoptic map is pasted together from 27 days 
of data (because the whole solar surface cannot be seen from 
Earth at a time). Additionally, the observational technique 
misses very thin structures. Keeping this in mind the match 
between the data sets seems reasonable. 
The location of the gridpoints of the map of magnetic field 
strength at the solar surface do not necessarily coincide with 
the derived footpoint locations and interpolation is applied, 
allowing to compute the magnetic field strength at the foot- 
points of the open field lines B(Rq, §o, ipo) and from this the 
expansion factor f s at every grid point on the source surface. 
The footpoint locations also allow to compute d F p for every 
corresponding grid point on the source surface by finding the 
nearest point that can be designated as part of a coronal hole 
boundary (yellow dots in Figure 3) . Tuning parameters have 
been adopted from the original paper (v m i n = 154 km/s, 
Vdei = 300 km/s, V exp = 0.3, a = 7.4 km/s, and /3 = 3.5). 
The resulting initial radial velocities at the inner radial grid 
boundary are shown in the respective top left panels of Fig- 
ures 4 and 6 for both minimum and maximum conditions. 
Typical solar minimum conditions are rather well known 
(high-speed streams in polar regions and low-speed streams 
in the ecliptic) and this general structure is reconstructed 
here. An interesting feature of high-speed wind in the equa- 
torial region is also found, and a comparison with OMNIweb 
http://omniweb.gsfc.nasa.gov/ data at the time in ques- 
tion does indeed show large equatorial speeds at the Earth 
at the time in question (see Figure 5 showing wind speeds of 
up to 700 km/s), which supports the assumption that this 
feature is real. For solar maximum there is neither a typical 
structure nor a sufficient amount of out-of-ecliptic measure- 
ments to compare with. There are however attempts to 
describe solar maximum conditions [Zurbuchen et al., 2004], 
but a comparison is still difficult, because the background 
field is constantly distorted by the large number (oc 5/day) 
of CMEs (as is also pointed out in Zurbuchen et al. [2004]). 
We therefore omit a detailed comparison for solar maximum 
because CMEs are not incorporated in this model as of yet. 
Equation (20) uses the value of the radial velocity and in its 
original form assumes constant mass flux F maS s, thus yield- 
ing an inverse relation between mass density p and radial 
velocity v r . We calculated F ma ss ,o by taking values of pv r 
at the Earth at 1 AU and applying an r 2 scaling to get val- 
ues at 0.1 AU. For this, 27-day averages of both quantities 
at the time in question are taken from the OMNIweb inter- 
face. The adopted values can be found in Table 2. 

We introduced an improvement to be in agreement with 
Ulysses data [McComas et al, 2000], which showed for solar 



Table 2. 27-day averaged values for mass density p and ra- 
dial velocity v r at 1 AU taken from the OMNIweb interface. 



Figure 5. OMNI solar wind speed data around 1987.2 
measured with the IMP 8 spacecraft. 
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Figure 6. Same as Figure 4 for solar maximum conditions. 

minimum conditions that the mass flux associated with high- 
speed streams is about half the one of low-speed streams. 
This is implemented by having the mass flux depend on the 
radial velocity by introducing a transition function f(v r ) 
(see Scherer et al. [2010]) defined as 

/(«,.) = i(tanh(A(u r + v t )) - tanh(A(u r - v t ))) (28) 



00 




that gives a sharp, but steady transition at v t with A = 20. 
The value of v t is chosen such that the transition from high 
to low speed occurs at i? « ±35° for our solar minimum data, 
which yields Vt = 1.5vq. Mass flux can then be calculated 
via 

Fmass = (0.5 + 0.5f(v r ))F mass ,o- (29) 
Results are shown in Figures 4 and 6 (b). 
Equation (21) is the pressure balance between gas pressure 
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and magnetic pressure solved for proton temperature (which 
is the only temperature in this one- fluid model). The tune- 
able parameter p to t is the total pressure, but its optimal 
values derived by Detman et al. [2006] change over the solar 
cycle, so that a rough tuning has been performed on our 
part to get realistic results and we chose ptot,min = 3.8po 
and ptot,max = 4.2p for solar minimum and maximum, re- 
spectively. 

The radial magnetic field strength has to be scaled from 
source surface values B 3S (R sa ) to values at the lower ra- 
dial grid boundary at R g b. The scale factor simply mimics 
the usual r~ 2 behavior, so that b aca ie = (Rss / R g b) 2 ■ This 
is a very simple approach, which is not optimal everywhere, 
because an inhomogeneous azimuthal magnetic field compo- 
nent yields deviations from the r~ 2 behavior for the radial 
magnetic field because of the solenoidality condition. How- 
ever, at small radial distances, such deviations should not 
be too serious since azimuthal components are small. 

Equations (23) and (24) set the $ components of both the 
velocity and the magnetic field equal to zero. We found that 
both v$ and B$ remain negligible throughout the computa- 
tional volume, so that these approximations are affirmed. 

As was demonstrated in the WD test case, the Parker 
spiral can be realized in two frames of reference (i.e. one 
co-rotating with the Sun and the other at rest with respect 
to the Sun's rotation), both of which have been found to 
give equivalent results. Equations (25) and (26) model the 
Parker spiral in the co-rotating frame, whereas in Detman 
et al. [2006] the rest frame is implemented. The equation 
for the azimuthal velocity mimics the behavior of the low 
corona rotating approximately like a rigid body, thus using 
the solar rotation frequency Q multiplied with the radius of 
co-rotation, which Detman et al. [2006] put at R v = 1.5.R®. 
This approach is questionable in view of the behavior of 
the azimuthal velocity derived by Weber and Davis [1967] 
shown in Figure 1, because from the radius of co-rotation 
outwards an approximate r~ x behavior is evident that would 
yield considerably smaller values at R g b- Correction for this 
with a respective scale factor has been applied here, so that 
R^ = 1.5-R©/(7? 9 f,/ 'Rq). This has been done in the recent 
update in Detman et al. [2011] as well, as R v is now de- 
scribed as "a length scale which determines the azimuthal 
velocity", but no specific value is given to compare with. 
This is still a rather crude estimate since the radius of co- 
rotation depends on the magnetic field strength near the Sun 
and simple 1/r scaling does not exactly reproduce the az- 
imuthal velocity solution that should be similar to the WD 
solution. But for this simple empirical model it shall suffice. 
The equation for the azimuthal magnetic field (26) is the 
result of the induction equation in the rotating steady state 
corona in corotating spherical coordinates. This is a dif- 
ferent approach then the one used for the testcases (see 
Equation (17)), which became necessary for the following 
reason: In the Weber-Davis-model the azimuthal magnetic 
field is homogeneous and can be extrapolated into the inner 
boundaries ghost cells from the values obtained from the 
innermost computed cells. The observationally based input 
for the magnetic field used here, however, is inherently in- 
homogeneous and makes extrapolation a challenging task, 
because of the staggered mesh implemented in CRONOS. 
Occasionally, this even gave rise to instabilities. To avoid 
this, we chose to fix the azimuthal magnetic field compo- 
nents according to Equation (26). In view of the inherently 
empirical nature of the whole set of equations to determine 
the remaining MHD quantities, this approach seems justifi- 
able. 

4.3. Radial initialisation 

The values of the MHD quantities at the lower radial 
boundary derived above will now also serve to find initial 
values for the rest of the numerical domain. Here, radial 



velocity along with mass density is to fulfill the stationary 
equation of continuity, which is achieved by having v r oc r 
and p oc r~ 3 . Azimuthal velocity in the co-rotating frame of 
reference is basically proportinal to radial distance (v v oc r). 
For adiabatic expansion the radial dependence of tempera- 
ture should approximately behave like T oc r~ 4 ' f where / is 
the number of degrees of freedom, related to the polytropic 
index via / = 2/(7 — 1). In our case of 7 = 1.05 this gives 
/ = 40. 

The magnetic field has to be initialized in a divergence-free 
manner, so that for later times the code retains a divergence- 
free field due to the constrained transport algorithm. This 
has been achieved by going through every cell layer radi- 
ally, and solving the solenoidality condition to find the next 
cells magnetic field components. A detailed description of 
this algorithm would require a discussion of the grid layout 
and the operation of the code, and is therefore omitted here 
[, but see Kissmann et al., 2009; Kissmann and Pomoell, 
2012]. 

4.4. Boundary conditions & Grid layout 

Boundary conditions are chosen similarly to those used 
for the test cases, i.e., reflecting 1} boundaries, periodic 
boundary conditions for 95, and an outflow condition at the 
outer radial boundary. At the inner radial boundary specific 
conditions are applied for each on quantity. 
The radial velocity and mass density were kept fixed at the 
inner radial boundary. Additionally, the ghost cells are as- 
signed values derived from extrapolation according to the 
governing power law in the first computed cell. Magnetic 
field components and v$ remain fixed to their initial values 
as well, but no extrapolation for the ghost cells is applied. 
The azimuthal velocity v v is expected to be dominated by 
the transformation term to the co-rotating frame of refer- 
ence, giving v v ~ — f2rsin(#). Thus, a radial extrapolation 
for the inner radial ghost cells is applied according to 

v Vli = — v v fi (30) 
ro 

where i indicates the negative index of the ghost cells, while 
the first cell in the computational domain carries the index 
0. 

For the observationally based input data the following grid 
layout has been chosen: For the latitudinal coordinate the 
ghost cells at each boundary should not coincide with or go 
beyond the poles since in the current setup coordinate sin- 
gularities have to be excluded, but this can be treated in the 
near future as well. Here, we are not primarily concerned 
with the polar regions. 

Furthermore, using the highest possible angular resolutions 
(1° x 1°) with an adequate radial resolution was found to 
require too long to converge. Instead we took half the pos- 
sible angular resolution and a comparison with full resolu- 
tion simulations yielded no significant differences. This gives 
tfgnd £ [8°, 172°] and a cell size of Ai? = 2° gives 172 cells. 
Similarly, for the longitudinal coordinate <p gri d £ [0°,360°], 
requiring 180 cells. The radial length of the simulation box 
is determined by Earth's orbital distance (215-Rq) and the 
lower radial grid boundary (at 21.5-Rq), where we chose a 
radial cell size of 1.5Rq. 



WIENGARTEN ET AL.: INNER-HELIOSPHERIC MHD SIMULATIONS 



X - 11 



(a) (b) 




Figure 7. Visualization of magnetic field topologies from the start surface out to 1 AU at solar minimum 
(a) and maximum (b): The spherical inner radial grid boundary is color-coded with the radial magnetic 
field B r , and the current sheet is the contour of B r — 0. Field lines are shown in black and demonstrate 
the bending caused by solar rotation. 



5. Simulation results 

5.1. Global magnetic field structure 

To get an impression of the results, pseudo three- 
dimensional visualizations of the respective data sets are 
shown in Figure 7. The inner radial grid boundary can be 
extrapolated from the set of data by performing a spher- 
ical clip at the respective radius (giving the sphere in the 
center). After applying 1/r 2 scaling, the colorbar is for the 
radial magnetic field strength B r , and is in accordance with 
the input data as was shown in Figure 2. The current sheet is 
the B r = contour through the computational grid. It can 
be seen that, at solar minimum, the current sheet is slightly 
wavy as expected from the input data, and it is subject to 
solar rotation as well, which can be seen by the bent wavy 
features. Magnetic field lines are colored black and are not 
uniformly distributed. The spiral structure can be made out 
as well by the bent field lines. For solar maximum condi- 
tions the current sheet's structure is rather extreme, with 
steep gradients in excursions towards the poles. 



5.2. Plasma structure at 1 AU 

The results at the outer radial grid boundary for radial 
velocity, mass density and magnetic field components are 
shown on the right hand sides of Figures 4 and 6, respec- 
tively. 

Comparing with the initial data, it can be seen that the 
results are strongly dependent on the input data, i.e. the 
input determines the results. Initial features can still be ob- 
served at the outer radial boundary but appear smoother. 
The effect of solar rotation is also visible as features shift 
in longitude. At solar minimum, for velocity and mass den- 
sity it can be observed that the transit from slow to fast 
speed is sharper than it was initially due to the differently 
strong acceleration. In equatorial regions we get speeds of 
about 400 km/s and the high-speed feature is at 700 km/s 
in accordance with the observational data taken from the 
OMNIweb (Figure 5). At polar regions, speeds of about 
650 km/s are found. Since there is no observational data for 
polar regions from that time it cannot be judged whether 



this value agrees with the conditions present at that time. 
Mass density is inversely correlated with radial velocity and 
the obtained values at Earth's orbit again agree with the 
observational data. The features of the magnetic field input 
are also preserved but the current sheet slightly broadens, 
which may be due to the finite spatial resolution. An in- 
teresting feature that is not present initially occurs at the 
location of the high-speed feature where we find an enhance- 
ment of the magnetic field strength. This is associated with 
the compression of the magnetic field when fast solar wind 
runs into slow one ahead, also known as co-rotating interac- 
tion regions. 

At solar maximum there is no typical global structure in 
wind speed features to compare with so that with the lim- 
ited number of observational data it is difficult to determine 
the quality of our results. It can, however, be stated (two 
top-right panels of Figure 6) that the obtained values for 
wind speed and mass density are again in the right order 
compared to the observational data. The initial magnetic 
field is largely preserved but is again modulated by features 
in the wind speed. 

5.3. 3-D velocity structure 

The radial velocity profiles at solar minimum are shown in 
Figure 8 in a meridional (ip = 0) and an equatorial (i9 = tt/2) 
slice through the computational grid. The meridional slice 
illustrates the dependence on latitude and Figure 9 shows a 
line plot accompanying the meridional slice, in which the ra- 
dial velocity curves for different latitudes can be seen. The 
top curves correspond to high-speed polar regions and the 
bottom ones to equatorial regions. The red-dashed line is the 
Parker solution for a temperature T = 0.49T = 1.4 • 10 6 K 
corresponding to the initial temperatures at polar regions, 
and shows that our results are close to it. The drop in some 
of the curves is due to slow wind features shifting to this lon- 
gitude at larger radial distances, which is an effect of solar 
rotation. These effects are also visible in the right panel of 
Figure 8. At close inspection it can also be seen that the shift 
in longitude between the start surface and the outer radial 
grid boundary due to solar rotation depends on the radial 
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Figure 8. (a) Meridional (ip = 0) and (b) equatorial slice through the computational grid for radial velocity. 



velocity. This corresponds to the so-called Parker angle that 
is the ratio between azimuthal and radial components. At 
the Earth it is usually at about 45° for radial velocities of 
about 400 km/s, which is also the case here near the ecliptic 
and we, too, find the corresponding shift in longitude to be 
at about 45°. We compared the winding angle (i.e. the ra- 
tio between the magnetic field components) with the Parker 
prediction and found good agreement for all latitudes except 
at the current sheet, which is understandable, because the 
Parker prediction does not involve dipolar fields. 

6. Conclusions and outlook 

In this work we successfully used observationally based 
input data of the radial magnetic flux to drive an inner he- 
liospheric MHD simulation to calculate solar wind quantities 
at radial distances up to 1 AU. Testing the code has been 
done by a comparison to the simplified analytical models 
by Parker [1958] and Weber and Davis [1967]. It could be 
demonstrated that the code succesfully reproduced the an- 
alytic results. 

The input data for the radial magnetic flux and its deriva- 
tion are described in Jiang et al. [2010], which amongst sev- 
eral advancements allows for higher angular resolution than 
that of comparable observational data from synoptic mag- 
netograms. While here the focus lay on two representative 
sets of data, those being one set at solar minimum and one 
at solar maximum, data are available for a whole time series 
going back to 1976. The possibility to run the code for more 
sets of data is present, including the opportunity to combine 
them to drive a time-dependent simulation. 
The MHD code requires input data for all plasma quantities, 
and it is necessary to use empirical formulas to derive them. 
The respective set of equations was taken from Detman et al. 
[2006, 2011], who work on a very similar goal, i.e. space 
weather forecasting, for which they use input data similar 
to those used here. These kinds of interfaces are constantly 



Figure 9. Radial velocity profiles at different latitudes 
accompanying the meridional slice of Figure 8. The top 
curves correspond to high-speed polar regions and the bot- 
tom ones to equatorial regions. The red-dashed line is the 
Parker solution for a temperature T = 0.49Tb = 1.4-10 6 K 
corresponding to the initial temperatures at polar regions, 
and shows that our results are close to it. 



augmented as testing goes on. In this respect, this work can 
be understood to be part of this process with the difference 
of using more sophisticated input data and applying the em- 
pirical formulas developed and tuned for a specific code to 
the one used here. 

The results give values for solar wind quantities at 1 AU that 
can be used as input for outer-heliospheric large-scale mod- 
els [e.g. Ferreira et al., 2007; Scherer et al., 2008], as back- 
ground solar wind for coronal mass ejection simulations [e.g. 
Kleimann et al, 2009] or to model solor wind plasma struc- 
tures [e.g. Dalakishvili et al., 2011, and references therein] 
such as co-rotating interaction regions. 
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